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Status  of  Effort: 


The  objective  of  this  project  is  the  development  and  study  of  a  new  class  of  propulsion  device  call 
the  Microdischarge  Plasma  Thruster  (MPT).  The  MPT  thruster  is  intended  for  small  satellite  propulsion 
and  is  expected  to  have  several  advantages  over  competing  small  satellite  propulsion  devices,  such  as  cold 
gas  microthrusters,  MEMS  fabricated  resistojets,  microPPT,  FEEP  devices,  and  miniaturized  version  of 
Hall  thrusters  and  Ion  thrusters.  The  primary  advantages  of  the  MPT  being  operation  with  a  variety  of 
propellant  gases  (noble  and  molecular  gases),  very  small  dimensions,  and  low  voltage  (~  lOCTs  V)  direct- 
current  operation. 

The  3-year  effort  discussed  in  this  report  is  broadly  divided  into  an  experimental  component  and 
a  computational  modeling  component.  We  have  achieved  a  reasonable  degree  of  success  with  both 
components  of  our  research  during  this  project. 

As  will  be  discussed  in  the  experimental  section,  this  project  has  enabled  us  to  conclusively 
demonstrate  that  it  is  possible  to  establish  a  stable  and  self-sustained  microdischarge  plasma  in  a  gas 
propellant  flow-through  situation  with  gas  efflux  into  vacuum.  A  number  of  diagnostics  have  been 
implemented  including  an  ion  beam  probe,  and  emission  spectroscopy  to  determine  the  physical  and 
chemical  state  of  the  microplasma  and  the  plume  from  the  MPT  device.  We  have  also  identified  and 
discussed  the  problem  of  erosion  in  the  device.  Details  of  the  experimental  are  given  in  the  next  section. 

In  the  computational  work  we  have  made  significant  progress  in  developing  an  entire  suite  of 
tools  to  provide  high-fidelity  computations  of  the  MPT  thruster.  A  self-consistent,  multi-species,  multi¬ 
temperature  model  for  the  MPT  discharge  and  a  model  for  the  low-Reynolds  number,  high-Mach  number 
compressible  flow  through  the  MPT  have  been  developed  and  coupled  together.  The  models  have  been 
exercised  to  provide  a  detailed  understanding  of  the  MPT  propulsion  device  operation.  Details  are 
provided  in  the  next  section. 
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Accomplishments  and  New  Findings: 

Experimental  Investigation  of  Micro  Plasma  Thruster  (MPT) 


I.  Introduction 

In  comparison  to  large  satellites,  small  satellites  are  more  suitable  choice  for  formation  flying  and 
multivehicle  operations  [1].  Because  of  small  size  and  mass,  they  are  also  relatively  easier  to  test  and 
assemble.  Such  satellites  typically  weigh  less  than  100  kg.  In  this  category,  nanosatellites  weigh  -10  kg 
and  have  linear  dimensions  of  less  than  50  cm.  Their  thrust  requirements  range  from  -lpN  for  small 
impulse  bit  precision  applications  to  -10  mN  for  satellite  slew  maneuvers.  The  construction  of 
nanosatellite  entails  miniaturizing  its  components  like  the  power  source,  electronics,  navigation  system, 
propulsion  system,  etc.  However,  a  scaled  down  version  of  a  large  pow  er  source,  electronics  or  a 
propulsion  system  is  not  always  viable,  and  so  alternative  approaches  must  be  explored  to  invent  smaller 
systems. 

There  is  significant  ongoing  research  in  the  space  propulsion  community  to  develop  micropropulsion 
devices,  but  options  for  nanosatellites  are  still  very  limited.  Among  the  electric  thrusters,  the  miniature 
hall  thruster  [2],  [3],  ion  thruster  [4],  [5]  and  micro  pulsed  power  thruster  (microPPT)  [6]  are  some 
examples  of  scaled  down  versions  of  large  thrusters.  The  miniaturized  hall  thruster  and  ion  thruster 
require  operating  power  of -100  W,  which  is  in  excess  of  what  can  be  accommodated  in  a  nanosatellite 
system  where  bus  powers  are  less  than  10-50  W.  In  contrast  to  this,  the  micro  pulsed  power  thrusters 
(microPPT)  can  operate  at  -10  W,  but  a  peak  pulse  voltage  can  be  several  kV  which  imposes  severe 
constraints  on  reducing  the  size  of  the  power  processing  unit  (PPU).  Thrusters  that  are  specifically  geared 
towards  micro  propulsion  also  suffer  from  high  voltage  requirement.  For  example,  the  field  emission 
electric  propulsion  (FEEP)  thruster  has  high  specific  impulse  (-10,000  s)  and  it  can  generate  thrust  levels 
of  (10  pN),  but  the  input  voltage  is  -10  kV  [7].  An  additional  drawback  of  the  FEEP  thruster  is  that  it 
uses  a  liquid  metal  propellant  like  the  liquid  Cesium,  which  can  pose  contamination  hazard.  Colloid 
thruster  is  also  similar  to  FEEP  thruster  [8].  Its  performance  characteristics  are  comparable  to  that  of  the 
FEEP  thruster;  specific  impulse  is  -1000  s  and  thrust  level  is  -  100  pN,  but  it  also  has  high  voltage 
requirements.  In  this  way,  although  some  electric  micro  thruster  devices  have  been  proposed  they  do  not 
provide  complete  solution  for  nanosatellite  propulsion,  mainly  because  of  their  high  operating  voltage  and 
power. 

Small  scale  versions  of  cold  gas  thrusters  [9],  resistojets  [10]  and  combustion  based  thrusters  [11]  are 
also  not  desirable  for  nanosatellite  propulsion.  In  combustion  based  micron-scale  devices  it  is  difficult  to 
stabilize  combustion  and  store  reactive  propellant.  The  issues  related  with  preheating  and  lifetime  of 
catalysts  can  also  complicate  combustion  in  a  micron-scale  devices.  The  miniaturized  cold  gas  thrusters 
and  resistojets  have  low  specific  impulse  (-  100  s  or  less),  which  is  also  not  ideal  for  micropropulsion. 

High  specific  impulse  can  possibly  be  attained  in  a  scaled  down  version  of  plasma  based  arcjet 
thruster,  but  miniaturization  of  such  a  thruster  is  limited  by  plasma  quenching  and  required  conditions  for 
plasma  breakdown.  As  predicted  by  Paschen  curves  [12],  the  minimum  value  of  breakdown  voltage  for 
most  gases  is  -100  V,  and  at  the  minimum  break  down  voltage,/?^,  which  is  the  product  of  pressure  (p ) 
and  inter  electrode  distance  ( d ),  is  -10  Torr-cm.  The  breakdown  voltage  increases  sharply  for  higher  or 
lower pd  values.  As  the  discharge  dimension  or  d  decreases  the  required  break  down  voltage  can  become 
extremely  high.  Also,  when  the  size  of  the  propulsion  device  is  decreased  the  ratio  of  its  surface  area  to 
volume  increases.  For  a  given  volume  of  plasma,  this  provides  larger  area  for  charged  species  to 
recombine,  and  so  the  quenching  phenomenon  hinders  plasma  ignition  and  stabilization  in  a  micron-scale 
device. 

The  MPT  overcomes  plasma  quenching  and  high  voltage  requirements.  Unlike  an  arc  jet  thruster,  it  is 
based  on  a  direct  current  micro  discharge.  The  micro  hall  thruster  (MHT)  [13]  also  utilizes  a  direct  current 
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micro-discharge,  but  it  requires  an  external  magnetic  field  to  stabilize  the  plasma.  The  plasma  in  an  MPT 
is  self  sustaining  and  so  an  external  magnetic  field  is  not  necessary. 

As  shown  in  Fig.  1,  the  MPT  has  a  multilayer  sandwich  structure.  It  has  three  electrodes-  eb  e2,  e3 
which  are  -  100  pm  thick  and  they  are  separated  by  -200  pm  thick  dielectrics-  m^  m2  and  m3.  As 
illustrated  in  Fig.  1  it  has  a  narrow  duct  -100  pm  in  diameter  which  expands  to  -330  pm  downstream  of 
the  throat.  From  gas  dynamic  point  of  view  the  cavity  inside  the  MPT  resembles  a  converging  diverging 
nozzle  that  can  accelerate  the  gas  to  super  sonic  velocities  downstream  of  the  throat. 

The  operating  mechanism  of  the  MPT  is  as  follows.  The  discharge  is  generated  between  ej  and  e3.  It 
heats  the  neutral  gas  which  accelerates  as  it  flows  in  the  converging  diverging  nozzle.  The  ion  flow  is  also 
accelerated  by  the  nozzle  and  the  floating  potential  that  exists  between  the  MPT  and  the  vacuum.  The 
expulsion  of  hot  gas  and  ions  into  the  vacuum  generates  thrust.  So  far  the  MPT  has  been  operated  using 
argon.  When  the  argon  flow  rate  is  5  seem,  the  pressure  in  the  cavity  ranges  from  -100  Torr  upstream  to 
-1  Torr  near  the  exit.  Between  e2  and  e3,  because  the  number  density  of  gas  species  (TV)  is  low,  the  ratio  of 
characteristic  electric  field  strength  (E)  to  the  number  density  of  neutral  species  EJN  is  high.  This  is 
favorable  to  enhance  ionization  process  in  the  discharge.  Ordinarily,  at  pressures  as  low  as  -1  Torr  and 
linear  dimensions  of  the  order  of  few  hundred  microns,  the  break  down  voltage  is  in  the  order  of  several 
kV  if  the  discharge  is  to  be  sustained  only  by  secondary  electron  emissions  at  the  cathode.  The  structure 
of  the  MPT  helps  to  overcome  such  high  voltage  requirement  to  sustain  the  discharge.  The  high  pressure 
(-100  Torr)  betw  een  ei  and  e2  helps  to  bring  down  the  break  down  voltage  of  gas  inside  the  MPT.  The 
region  of  the  discharge  between  ei  and  e2  can  also  supply  electrons  to  the  region  between  e2  and  e3  to 
increase  and  sustain  ionization  process  in  that  region.  In  this  way,  MPT  can  be  stably  operated  using  a 
non  corrosive/  non  contaminating  propellant  at  relatively  low  voltages. 


II.  Experimental  Set  Up 


The  MPT  was  operated  inside  a  vacuum  chamber.  Experiments  were  carried  out  to  determine  the 
stable  operating  conditions,  electrical  characteristics  of  MPT  and  properties  of  the  plume.  The  following 
section  describes  the  apparatus  and  experimental  set  ups. 

A.  Vacuum  System 

The  MPT  is  operated  in  a  non-magnetic  stainless  steel  vacuum  chamber  which  is  0.15  m  in 
diameter.  The  chamber  consists  of  six  symmetrically  placed  extensions  of  length  0.08  m  and  diameter  0.1 
m,  which  provide  access  to  the  interior  of  the  chamber.  Optical  access  is  provided  by  the  glass  viewports 
on  the  extensions.  The  glass  is  4  mm  thick  and  of  type  kodial,  which  means  high  frequencies  like  UV  are 
absorbed  by  these  viewports.  The  chamber  is  evacuated  using  a  Varian  TV-551  turbo  molecular  pump 
backed  by  Alcatel  2015  SD  general  vane  roughing  pump.  Without  gas  flow,  when  the  roughing  pump 
evacuates  to  15  mTorr,  the  turbo  molecular  pump  is  switched  on  to  lower  the  chamber  pressure  to  -  1 
pTorr.  WTen  the  pressure  in  the  chamber  is  between  atmospheric  and  15  mTorr  it  is  measured  by 
KJL912161  convection  gauge,  which  has  a  measurement  range  of  lxlO'3  Torr  to  760  Torr.  When  the 
turbo  molecular  pump  also  starts  to  operate  the  chamber  pressure  goes  below  1  mTorr  and  so  it  is 
measured  using  an  FIPS  943  digital  Penning  ionization  cold  cathode  gauge,  which  has  measurement  range 
of  lxlO'10  Torr  to  lxlO'2  Torr  .  The  cold  cathode  gauge  is  mounted  on  the  top  extension  of  the  chamber. 
The  distance  between  this  pressure  gauge  and  the  center  of  the  chamber  is  0.25  m.  The  reservoir  pressure 
is  measured  0.62  m  upstream  of  the  chamber  by  an  MKS  Baratron  Capacitance  manometer.  The 
propellant  gas  flow  rate  is  controlled  by  an  Omega  FMA  2615 A  mass  flow  controller.  The  maximum 
possible  flow  rate  is  5.120  seem.  Without  gas  flow  the  back  pressure  is  10"6  Torr,  and  the  reservoir 
pressure  is  0.9  Torr.  When  the  gas  flow  rate  is  1-5  seem  (17-83  pL/s),  the  reservoir  pressure  is  10-100 
Torr,  and  the  back  pressure  is  of  the  order  of  lO'4  Torr. 
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B.  MPT  Device 


The  schematic  of  MPT  is  given  in  fig.  1.  It  is  fabricated  by  the  experimenters.  The  electrodes  are 
cut  out  from  a  sheet  of  99.95  %  Molybdenum  foil  that  is  100  pm  thick.  The  upstream  electrode  (ei)  and 
the  middle  electrode  (e2)  are  25.4  mm  x  6  mm.  The  downstream  electrode  (e3)  is  25.4  mm  x  20  mm.  The 
dielectric  is  cut  out  from  mica  sheet  that  is  200  pm  thick.  The  dielectric  layers  are  25.4  mm  x  25.4  mm. 
The  diameter  of  the  cylindrical  cavity  is  100  pm  from  mi  to  m3,  and  it  is  330  pm  in  e3.  The  flute  length  of 
the  100  pm  drill  is  about  2  mm  and  it  can  break  easily  when  tried  to  penetrate  more  than  1  mm.  The 
drilling  is  carried  out  at  rpm  of  2800.  Higher  rpm  is  desirable  as  the  drill  bit  is  thin.  First  the  layers  from 
m:  to  m3  are  assembled.  They  are  held  together  by  epoxy  resin.  Then  a  bore  is  drilled  through  this  piece 
using  a  drill  bit  of  diameter  100  pm.  The  layer  e3  is  drilled  separately  using  330  pm  drill  bit.  Layer  e3  and 
m3  are  attached  using  RTV  108.  The  holes  in  e3  and  m3  are  aligned  by  viewing  them  under  a  x10 
magnifier.  The  mass  of  the  assembled  piece  is  0.9  g  and  the  total  thickness  is  about  1.1  mm. 

As  shown  in  fig.  2,  after  machine  work  the  MPT  is  attached  on  a  lexan  block  using  screws  and 
nuts.  A  tube  of  6  mm  diameter  feeds  argon  to  the  cavity  inside  the  MPT  through  a  hole  in  the  lexan  block. 
The  hole  in  the  lexan  block  that  overlaps  with  100  pm  hole  in  mi  is  0.980  mm  in  diameter.  The  exposed 
surfaces  of  the  electrodes  and  wires  are  insulated  with  kapton  tape  to  prevent  unwanted  arcing. 

C.  Electrical  Circuit 

The  schematics  of  circuits  used  to  operate  the  MPT  are  shown  in  Fig.  3.  Stable  operation  is 
possible  in  two  possible  ways.  As  shown  in  fig.  3  a,  ey  is  the  anode,  e3  is  the  cathode  and  e2  is  at  a 
potential  between  e*  and  e3.  In  fig.  3  b,  ei  is  anode  and  e3  is  cathode  whereas  e2  is  floating.  The  power 
supply  for  ei  is  Matsusada  AU-2P1 50-  LC.  It  can  supply  maximum  voltage  of  2  kV  and  maximum  current 
of  150  mA.  The  power  supply  for  e2  is  Matsusada  AU-3R-100.  It  can  also  supply  maximum  voltage  of  2 
kV,  but  minimum  current  of  100  mA.  Both  power  supplies  can  be  set  to  generate  positive  or  negative 
voltages.  The  values  of  ballast  resistors,  RBi  is  160  kO  and  RB2  is  60  kf2.  The  resistances  of  current  view 
resistors  (CVR)  are  lQkeach.  Voltages  across  CVR’s  are  measured  using  Radio  Sha  ck  digital 
multimeters. 

D.  Ion  Probe  Set  up 

An  ion  probe  was  used  to  determine  ion  distribution  in  the  plume.  The  set  up  is  shown  in  fig.  4. 
As  shown,  the  probe  consists  of  a  copper  rod  of  radius  1  mm  inserted  into  a  glass  tube  that  has  1  mm 
thick  wall.  The  probe  surface  is  biased  at  -100  V  using  Matsusada  AU-3R100.  As  shown  in  fig.  5,  the 
fringing  effect  is  negligible  at  a  bias  voltage  of  -100  V.  This  is  determined  by  checking  the  variation  in 
ion  current  as  the  bias  voltage  is  made  increasingly  negative.  As  the  probe  is  made  more  negative  it  also 
pulls  ions  that  may  otherwise  be  incident  outside  its  periphery,  and  so  although  the  plume  current  is 
constant  the  probe  records  higher  current  when  its  negative  bias  is  increased. 

The  probe  is  mounted  on  a  linear  translation  feed  through  so  that  it  can  traverse  perpendicular  to 
the  plume  axis  on  the  horizontal  plane.  When  the  probe  aligns  with  the  plume  axis  it  is  at  4.4  cm  from  the 
MPT  exit.  A  1  kQ  CVR  and  a  Radio  Shack  digital  multimeter  were  used  to  read  the  ion  current. 

E.  Optical  Imaging  Setup 

Images  of  the  plume  were  captured  using  a  Cannon  PowerShot  A95  camera.  Typical  ISO  value  was 
200,  and  F  number  was  3.2.  The  camera  was  positioned  outside  the  chamber  at  about  0.4  m  from  the 
plume.  Power  of  the  MPT  and  the  exposure  time  of  the  camera  were  varied  during  the  experiment.  After 
an  MPT  burnt  out,  layers  e3  and  m3  were  viewed  under  a  CZM4  microscope  to  study  the  damaging  effects 
of  the  hot  plume  on  the  MPT. 

F.  Spectroscopy  Setup 
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Optical  Emission  Spectroscopy  was  used  to  determine  the  electronic  and  rotational  temperatures  of  the 
plume.  Emission  from  the  plume  was  collected  outside  the  chamber  by  a  1 9  cm  focal  length  lens  which 
was  positioned  about  38.7  cm  from  the  exit  of  MPT.  The  diameter  of  the  lens  was  6.7  cm.  As  shown  in 
Fig.  4,  the  light  from  the  lens  was  focused  on  an  optical  fiber.  It  is  a  Thorlabs  M29L01  fiber,  which  has 
numerical  aperture  (NA)  of  0.39  and  core  diameter  of  600  jam.  The  light  from  the  other  end  of  the  optical 
fiber  was  collected  by  a  lens  which  had  focal  length  of  1.5  cm  and  diameter  of  2.8  cm.  It  was  then  focused 
on  the  entrance  slit  of  a  Newport  74050  spectrometer,  which  had  a  grating  of  1800  g/mm.  The  slit  width 
was  set  to  200  jam.  The  detector  mounted  on  the  spectrometer  is  a  Newport  78105  array  detector,  which  is 
a  3000  element  linear  CCD  detector  with  pixel  size  200  jam  (H)  x  7  pm  (W). 

The  depth  and  width  of  field  were  determined  by  probing  emission  from  an  argon  calibration  lamp 
focused  into  a  pinhole  that  had  a  diameter  of  100  pm,  and  they  were  found  to  be  3.03  cm  and  360  pm 
respectively.  The  distance  between  the  source  and  the  lens  was  approximately  41  cm.  The  lens  and  optical 
fiber  were  mounted  on  a  translation  stage,  which  could  move  them  to  the  left  and  right  (x-axis)  of  the  pin 
hole  or  towards  and  away  (y-axis)  from  the  pin  hole.  The  relative  distance  between  the  lens  and  the  fiber 
is  fixed  as  they  are  attached  to  the  translation  stage  so  that  it  displaces  the  lens  and  the  fiber  in  unison 
relative  to  the  source.  The  least  count  of  the  translation  stage  was  20  pm.  The  data  obtained  from  this 
procedure  are  shown  in  figs.  7. 

Fig.  7a  shows  the  counts  recorded  from  the  pin  hole  as  the  translation  stage  carries  the  lens  and  fiber 
along  the  x-axis.  ‘Integrated  counts’  is  calculated  by  integrating  the  area  under  the  line  emission  at  81 1 
nm.  Gaussian  fit  is  applied  to  the  data  shown  in  fig.  7a  to  estimate  the  FWHM.  The  probe  size,  which  is 
the  resolution  of  the  width  of  region  viewed  by  the  optical  system,  is  defined  to  be  equal  to  the  FWHM  of 
the  curve  fit.  Hence,  the  probe  size  is  360  pm. 

Fig.  7  b  shows  the  counts  recorded  from  the  pin  hole  as  the  collection  optics  is  moved  in  the  y- 
direction  from  the  pin  hole.  Lorentzian  fit  is  performed  on  the  data  in  fig  7  b  to  estimate  the  FWHM, 
which  is  defined  as  the  depth  of  field  of  the  optical  system.  From  the  curve  fit,  the  depth  of  field  is  3.03 
cm. 

The  plume  spectra  were  recorded  from  the  region  2  mm  downstream  of  the  exit  plane.  There  the 
length  of  the  visible  plume  perpendicular  to  its  axis  is  estimated  to  be  less  than  3.03  cm.  This  means  the 
lens  collects  light  over  the  entire  depth  of  the  plume  during  spectroscopy  experiments. 

The  recorded  intensities  of  the  emissions  were  calibrated  using  a  Cl-system  SR-20  cavity  blackbody 
that  can  reach  maximum  temperature  of  1200  °C.  The  wavelengths  of  the  spectra  were  calibrated  by 
recording  known  lines  from  Argon  and  Mercury  lamps.  As  shown  in  fig.  8,  according  to  546  nm  line 
recorded  from  the  Mercury  lamp,  the  spectral  resolution  is  0.452  nm. 


III.  Results  and  Discussion 

Stable  operation  modes  of  the  MPT  were  determined.  Its  current,  voltage  and  power  characteristics 
were  recorded.  Visual  features  of  the  plume  were  also  studied.  The  results  also  helped  to  understand  the 
factors  that  lead  to  deterioration  of  the  MPT  during  operation.  The  data  was  also  used  to  infer  the 
composition  of  the  plume  and  its  electronic  excitation  temperature.  The  following  section  presents  the 
results  and  their  interpretation. 


G.  Optical  Imaging 

As  shown  in  fig.  9,  the  plume  rapidly  expands  into  the  chamber.  The  pictures  are  taken  with  2.5  s 
exposure  time  when  the  argon  flow  rate  was  5  seem  for  input  powers  of  40  mW,  291  mW  and  531  mW. 
There  is  a  structure  in  the  appearance  of  the  plume.  It  has  a  luminous  core  anchored  at  the  exit.  The  core 
typically  extends  about  1.5  mm  downstream,  and  then  transitions  to  a  faint  bluish  region.  The  variation  in 
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the  brightness  of  the  plume  over  its  length  is  affected  by  its  expansion  and  the  lifetime  of  the  excited 
species  in  the  plume  In  addition,  the  presence  of  a  shock  wave  may  also  shape  its  visual  structure. 

The  MPT  is  operated  as  a  highly  under-expanded  jet  as  the  pressure  at  exit  is  ~  1  Torr  and  the 
pressure  in  the  chamber  is  -10''  Torr,  and  so  the  plume  is  expected  to  have  a  barrel  shock.  The  luminous 
core  of  the  plume  could  be  the  region  inside  the  barrel  shock.  The  transition  from  luminous  core  to  the 
faint  bluish  region  could  be  consequences  of  the  presence  of  the  Mach  disk  and  oblique  shocks.  The 
pressure  and  density  rise  behind  a  shock.  This  could  increase  the  rate  of  quenching  collisions.  This 
implies  that  the  excited  species  are  quenched  behind  the  shock  and  they  stop  emitting  which  leads  to  a 
faint  region  downstream  of  the  luminous  core. 

Schlieren  imaging  using  3  ft  focal  length  mirror  was  attempted  on  the  plume,  but  it  did  not  reveal 
the  expected  density  gradients.  Schlieren  imaging  could  not  even  distinguish  the  plume  from  the 
background  gas.  Clearly  the  set  up  that  was  used  was  not  sensitive  enough.  The  lack  of  sensitivity  of  the 
set  up  that  was  used  can  be  understood  from  the  formula  for  relative  intensity  variation.  The  relative 
intensity  variation  in  Schlieren  method  is  given  by 


AP_=A_  f  1  dn(x,y,z)  d_ 

p>  a  Ln  dy 

(1) 

,  where 

P,  =  Ifia 


,  and 

A P  =  IfiAa . 

I,  is  the  intensity  of  the  image  at  the  focal  plane  in  the  absence  of  plume,  and  b9  a  are  the  width  and  height 
of  the  image  at  the  focal  plane;  Aa  is  the  size  of  deflection  of  image  in  the  focal  plane  because  of 
refractive  index  of  the  plume;  f2  is  the  focal  length  of  the  mirror  that  focuses  light  on  the  knife  edge;  n  is 
the  refractive  index  of  the  plume  which  is  related  to  density  by 


N=l+Kp 

(2) 

Here,  K  is  called  the  Gladstone  constant  and  p  is  the  density  of  the  medium.  The  lack  of  sensitivity  of 
Schlieren  method  for  plume  diagnostics  arises  mainly  from  the  smallness  of  AZ ,  which  is  the  depth  of  the 
plume  along  the  optical  path.  The  sensitivity  will  also  improve  by  using  a  mirror  of  longer  focal  length. 

The  change  in  the  length  of  the  luminous  core  as  the  MPT  power  and  flow  rate  are  varied  is 
shown  in  fig.  10.  The  length  was  determined  by  multiplying  the  length  of  the  luminous  section  in  the 
photograph  by  the  ratio  of  the  real  length  of  the  lexan  block  to  its  length  in  the  photograph.  According  to 
this  data,  the  length  of  the  core  is  independent  of  the  flow  rate.  It  increases  almost  linearly  with  power. 

Images  of  e3  and  m3  were  taken  after  an  MPT  stopped  functioning.  As  shown  in  Fig.  11  they  reveal 
that  the  electrodes  melt  and  erode  during  operation.  The  melted  part  of  the  electrode  settles  on  the 
dielectric.  The  MPT  stops  to  function  when  the  meltdown  clogs  the  cavity.  The  shape  and  size  of  the 
holes  were  also  distorted  after  repeated  usage. 


H.  Operation  Modes 

The  performance  of  MPT  depends  on  the  relative  potentials  of  the  electrodes.  There  were  minimum 
fluctuations  in  the  electrode  currents  and  plume  intensity  when  ej  was  set  as  anode,  e3  as  cathode  and  e2 
was  floating.  When  anodes  and  cathodes  were  reversed  in  this  configuration-  that  is  e3  was  set  as  anode, 
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ei  was  set  as  cathode  and  e2  floating,  the  required  breakdown  voltage  was  about  1500  V;  when  the  plume 
appeared  it  quickly  damaged  the  MPT  and  made  it  defunct  so  that  it  was  not  possible  to  study  this 
configuration.  When  e2  was  also  powered  the  discharge  was  confined  between  e!  and  e2  and  a  plume  was 
not  observed.  After  repeated  operations  in  the  configuration  in  which  ei,  e2  and  e3  were  set  at 
monotonically  decreasing  potentials,  e2  and  e3  short  circuited  and  a  bright  plume,  which  was  relatively 
stable,  was  observed.  In  this  relatively  stable  configuration,  e2  was  effectively  floating  because  current 
flow  was  predominantly  in  ei  and  e3  only.  It  was  not  possible  to  generate  a  plume  by  an  external  short 
circuit  between  e2  and  e3.  We  believe  that  multiple  operation  of  MPT  melts  the  area  of  electrodes  exposed 
to  the  discharge  inside  the  MPT.  When  the  melted  part  settles  between  e2  and  e3,  the  e2  becomes 
effectively  floating;  current  flows  in  ei  and  e3  and  a  relatively  stable  plume  is  observed.  After  the  most 
stable  configuration  was  determined,  e2  was  deliberately  left  floating,  and  the  MPT  and  its  plume  were 
characterized.  In  this  configuration,  plume  can  be  generated  in  a  new  MPT  so  that  it  is  not  necessary  to 
operate  it  multiple  times  to  generate  the  plume.  In  other  configurations  it  was  not  possible  to  study  the 
plume  and  the  currents  in  the  electrodes  as  they  were  either  not  stable  or  the  break  down  voltage  was  too 
high. 

I.  Electrical  Characteristics 

When  ei,  e2  and  e3  are  set  at  monotonically  decreasing  potentials  their  currents  are  shown  in  fig. 
12.  In  this  configuration,  the  discharge  is  confined  between  ei  and  e2.  A  visible  plume  is  not  observed. 
Current  in  e3  (I3)  is  almost  zero.  Almost  all  current  flows  from  ei  to  e2.  When  the  voltage  is  further 
increased  the  plume  starts  to  flicker  and  I3  increases.  Currents  in  all  the  electrodes  (Ij,  I2,  I3)  become 
unstable,  and  it  is  no  longer  possible  to  record  the  current-voltage  characteristics.  Multiple  operations  lead 
to  an  internal  short  circuit  between  e2  and  e3.  fig.  13  shows  the  currents  when  e2  and  e3  have  an  internal 
short  circuit  and  e2  was  not  powered,  fig  14  shows  currents  in  ej  and  e3  when  e2  was  left  floating  and  there 
was  no  short  circuit  between  e2  and  e3. 

Ii  and  I3  are  higher  in  fig.  13  than  in  fig.  14  possibly  because  when  e2  and  e3  are  internally  short 
circuited  they  act  as  single  thicker  cathode.  The  currents  in  successive  runs  decrease  possibly  because  of 
the  erosion  of  cathode  during  operation.  When  all  three  electrodes  are  powered,  and  the  plasma  is 
confined  between  ei  and  e2,  the  erosion  of  e2  first  leads  to  internal  short  circuit  between  e2  and  e3.  This 
increases  the  current  as  the  effective  cathode  is  now  thicker.  After  current  starts  to  flow  in  e3,  its  erosion 
leads  to  decrease  in  current  values  in  successive  runs.  The  plots  show  that  the  typical  current  in  MPT  was 
-  0.1  mA  and  voltage  input  was  ~  1000  V. 

Fig.  15  shows  the  variation  in  the  power  of  MPT  as  supply  voltage  and  flow  rate  are  increased. 
Typical  operating  power  is  between  50  and  500  mW.  For  a  given  supply  voltage  the  power  increases 
when  the  flow  rate  increases.  When  flow  rate  increases  the  current  in  the  MPT  can  also  increase  because 
of  the  increase  in  the  flux  of  atoms  or  it  can  decrease  because  the  pressure  increases  with  flow  rate  and  an 
increase  in  pressure  increases  collision  frequency.  Current  density,  /  is  given  by 

J=  nq  Vd 

(3) 

=  nqjuE 

J=  nq2E/mv 

(4) 

,  where  n  is  the  number  density  of  charge  carriers,  q  is  the  charge,  Vd  is  the  drift  velocity,  ja  is  the 
mobility,  E  is  the  electric  field  strength,  m  is  the  mass  of  the  charge  carrier  and  v  is  the  momentum 
exchange  collision  frequency.  Hence,  J  is  inversely  proportional  to  v.  According  to  fig.  15,  in  MPT  an 
increase  in  flow  rate  at  a  given  supply  voltage  results  in  an  increase  in  current  and  power,  which  implies 
that  the  increase  in  the  number  density  of  charge  carriers  in  the  MPT  outweighs  the  increase  in  collision 
frequency. 

J.  Ion  distribution  in  the  plume 
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Measurements  carried  out  by  the  ion  probe  are  represented  in  fig.  16.  According  to  fig.  16  a,  the 
typical  ion  flux  is  IpA/cm  .  The  off-axis  angle  is  defined  as  the  angle  between  the  plume  axis  and  the 
probe’s  line  of  sight  of  the  MPT  exit.  Such  low  ion  content  in  the  plume  is  not  surprising  given  that  the 
power  input  is  250  mW.  It  can  be  shown  that  at  least  5.09  W  is  necessary  to  fully  ionize  Ar  flow  at  5 
seem  to  Ar+  even  if  all  the  input  energy  were  to  go  into  ionization. 

Fig.  16  b  contains  the  same  data  as  fig.  16  a  except  that  now  the  horizontal  axis  is  the  distance 
over  which  the  probe  travels.  Total  ion  current  in  the  plume  calculated  from  the  curve  fit  on  fig.  16  b 
is \1  .Q/uA  .  The  calculation  of  total  ion  current  assumes  that  the  ion  current  has  a  Gaussian  distribution 
over  the  cross  section.  It  can  be  shown  that  the  total  ion  current 

Assuming  that  all  the  ions  are  Ar+,  the  ratio  of  ion  flow  rate  to  flow  rate  of  neutral  species  in  the 
plume  is  5.2*  10‘5.  Assuming  that  all  the  ions  in  the  discharge  are  also  Ar+,  the  ratio  of  ion  current  in  the 
plume  to  discharge  current  turns  is  5.6 x  10'2. 

The  low  proportion  of  ion  content  implies  that  most  of  the  thrust  comes  from  gas  dynamic 
expansion  rather  than  ion  acceleration. 

K.  Spectroscopy 

Line  emissions  from  molybdenum,  argon  ions  and  excited  neutral  species  of  argon  were  detected 
between  405  nm  and  827  nm.  Mo  lines  were  predominant  between  500  nm  and  650  nm.  Ar+  lines  were 
identified  at  434  nm,  61 1  nm  and  617  nm.  Most  of  the  recorded  line  emissions  were  from  excited  neutral 
species  of  Ar  and  from  Mo.  The  presence  of  Mo  species  in  the  plume  implies  that  some  of  the  exposed 
areas  of  Mo  electrodes  are  evaporated  and  carried  away  by  the  plume.  It  was  not  possible  to  identify  all 
the  lines  that  were  recorded  because  the  spectral  resolution  was  0.452  nm.  Especially  at  wavelengths 
below  600  nm  there  are  relatively  strong  emissions  from  Mo  and  Ar+,  but  they  could  not  be  distinguished 
from  one  another.  Parts  of  spectrum  are  shown  in  fig.  3.12.  Doppler  fit  has  been  applied  to  some  of  the 
lines.  The  fitted  spectrum  is  used  to  calculate  the  intensity  from  a  given  line  emission  which  is  necessary 
to  calculate  the  electronic  excitation  temperature  as  explained  in  the  next  section. 

L.  Determination  of  temperature 

The  electronic  excitation  temperature  was  inferred  by  Boltzmann  fit  on  the  plot  of  Ln  [I/gAv]  vs. 
s/k,  where  I  is  the  intensity  of  the  line,  g  is  the  degeneracy  of  the  upper  level  of  energy  transition,  A  is  the 
Einstein  coefficient,  v  is  the  frequency  of  transition,  s  is  the  transition  energy  and  k  is  the  Boltzmann 
constant.  For  a  given  value  of  v  the  values  of  g,  A  and  s  were  obtained  from  the  NIST  database.  I  is 
calculated  by  integrating  the  area  under  the  Doppler  fit  on  spectral  lines.  The  gradient  of  line  fit  on 
Ln[I/gAv]  vs.  s/k  is  equal  to  -1/T,  where  T  is  the  electronic  excitation  temperature. 

Molybdenum  temperature  could  not  be  determined  as  the  energy  states  could  not  be  identified. 
The  energy  levels  of  metals  are  closely  spaced  and  it  is  often  difficult  to  distinguish  them. 

As  shown  in  figs.  18,  the  electron  temperature  inferred  from  data  on  neutral  species  is  1270  K, 
and  the  electron  temperature  inferred  from  ions  is  3510  K.  The  values  of  temperatures  obtained  here  are 
much  lower  compared  to  the  expected  electron  temperatures  inside  the  discharge,  where  the  charged 
species  are  accelerated  over  a  potential  difference  of -1000  V.  This  means  the  electronic  excitation  of  the 
atoms  and  ions  that  gave  rise  to  the  recorded  spectra  happened  in  the  after  glow  i.e.  outside  the  MPT.  As 
mentioned  above,  the  probe  region  was  located  2  mm  downstream  from  the  exit  plane.  Free  electrons 
cannot  reach  that  far  from  the  discharge  by  convection  so  the  inferred  temperature  cannot  possibly 
indicate  electron  population  distribution  at  the  probe  region.  The  ions  and  neutral  species  were  excited 
somewhere  upstream  of  the  probe  region,  and  they  reached  the  probe  region  by  convection  so  the  values 
of  temperature  given  here  provides  an  idea  of  electron  population  distribution  upstream  of  2  mm  in  the 
after  glow.  The  maximum  error  in  the  calculation  arises  from  the  noise  in  black  body  signal  used  to 
calibrate  the  intensity.  The  fact  that  the  line  fit  does  not  pass  between  all  the  maximum  and  minimum 
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values  of  /  A}ig}V  indicates  that  the  Ar  and  Ar+  species  in  the  plume  do  not  have  equilibrium 
distribution. 


IV.  Conclusion 

In  this  way,  MPT  is  an  electro-thermal  micro  propulsion  device,  which  was  designed  and  studied  with  the 
goals  of  minimizing  its  size,  mass,  operating  power  and  propellant  consumption.  The  MPT  was 
successfully  operated  at  low  flow  rates  of  ~  1  seem.  The  most  stable  operation  is  attained  when  the 
upstream  electrode  ei  is  anode,  the  downstream  electrode  e3  is  cathode  and  the  middle  electrode  e2  is 
floating.  The  MPT  operates  as  an  electrothermal  thruster  because  most  of  the  momentum  thrust  of  the 
MPT  is  generated  by  the  expulsion  of  hot  neutral  argon  atoms.  Measurements  from  ion  probe  revealed 
that  there  is  also  low  proportion  of  ions  in  the  plume;  therefore,  MPT  is  able  to  expel  ions  without  the  use 
of  an  external  magnetic  field  or  an  external  electron  emitter.  Comparison  of  current  voltage  characteristics 
of  the  MPT  with  the  ion  probe  measurements  show  that  the  net  ion  current  in  the  plume  is  also 
considerably  less  than  the  discharge  current.  The  ion  content  in  the  plume  was  also  detected  by  emission 
spectroscopy.  The  detection  of  Molybdenum  atoms  by  spectroscopy  implies  that  the  plume  melts  and 
evaporates  the  exposed  regions  of  the  electrodes.  The  magnified  view  of  the  exit  hole  in  e3  and  m3  also 
reveal  that  the  holes  are  distorted  and  enlarged  by  the  time  the  MPT  becomes  defunct.  The  electronic 
excitation  temperature  of  the  ions  and  neutral  species  were  found  to  be  3510  K  and  1270  K  respectively. 
The  thrust  from  the  MPT  is  yet  to  be  determined.  The  material  selection  issue  for  the  electrodes  also 
needs  to  be  addressed  in  future. 
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Electrode  1  Electrode  2  Electrode  3 


Figure  1.  The  schematic  of  MPT. 


(a)  (b) 


Figure  2.  (a)  Side  view  of  the  MPT  assembly  that  is  placed  inside  the  vacuum  chamber,  (b)  .  Front  view  of 

the  MPT  assembly. 
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Figure  3.  Electrical  circuit  of  MPT.  (a)  Middle  electrode,  e2  is  not  floating,  (b)  Middle  electrode  is  floating. 


Figure  4.  Ion  probe  set  up  inside  the  vacuum  chamber. 
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0 


Figure  5.  Fringing  effect  The  power  input  to  the  MPT  is  constant,  but  the  current  recorded  by  ion  probe 

increases  as  it  is  made  more  negatively  biased. 


Figure  6.  The  set  up  for  emission  spectroscopy. 
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Integrated  Counts 


(a)  (b) 

Figure  7.  (a)  Determination  of  probe  size.  Integrated  counts  is  the  area  under  the  line  emissions  at  811  nm 
wavelength  from  an  argon  calibration  lamp,  (b)  Determination  of  depth  of  field.  Integrated  counts  is  the 
area  under  the  line  emissions  at  811  nm  wavelength  from  an  argon  calibration  lamp. 


wavelength  (nm) 

Figure  8.  546  nm  line  recorded  from  mercury 
calibration  lamp.  FWHM  is  0.452  nm.  Slit  width  of  the 
spectrometer  was  200  pm. 
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(a)  (b)  (c) 

Figure  9.  (a)  Plume  image  with  2.5  s  exposure  time.  Input  power  is  40  mW,  and  the  luminous  length  is  0.56 
mm.  (b)  Plume  image  with  2.5  s  exposure  time.  Input  power  is  291  mW,  and  the  luminous  length  is  1.2  mm. 
(c)  Plume  image  with  2.5  s  exposure  time.  Input  power  is  531  mW,  and  the  luminous  length  is  2.2  mm. 
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Figure  10.  Variation  in  the  length  of  the  luminous  core  of 
the  plume  as  the  power  and  flow  rate  are  changed.  The 
length  of  the  core  increases  almost  linearly  with  the  power. 
It  is  independent  of  the  flow  rate. 
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currents  in  electrodes  (mA) 


Figure  11.  a)  m3  before  running  MPT.  b)  e3  after  MPT  stops  functioning,  c)  e3  before  running  MPT. 
d)  downstream  face  of  e3  after  MPT  stops  functioning,  e)  upstream  face  of  e3  after  MPT  stops 
functioning. 
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Figure  12.  MPT  current  vs.  voltage  when  all  three 

electrodes  are  powered.  Figure  13.  MPT  current  vs.  voltage  when  there  is  internal 

short  circuit  between  e2  and  e3. 
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Figure  15.  MPT  power  vs.  voltage  at  given  flow  rates. 

Figure  14.  MPT  current  vs.  voltage  when  e2  is  floating. 
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Figure  16.  (a)  Ion  flux  vs.  off-axis  angle. 
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(b)  Ion  flux  vs.  distance  travelled  by  the  probe. 
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Figure  17.  Some  spectra  from  plume  emissions.  The  Doppler  fit  is  used  to  calculate  the  total  intensity  from  a 

given  line  emission. 
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a)  b) 

Figure  18.  a)  Determination  of  electronic  excitation  temperature  of  Ar+.  b)  Determination  of  electronic 
excitation  temperature  of  Ar. 
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Computational  Investigation  of  Micro  Plasma  Thruster  (MPT) 


I.  Introduction 

Small  satellites  (mass  less  than  1 00  kg)  have  recently  gained  interest  for  various  commercial,  military, 
and  science  space  missions.1'1  Several  advantages,  including  a  reduced  launch  cost  are  realized  with 
small  satellites.  Small  satellites  are  power-limited,  typically  delivering  less  than  a  watt  per  kilogram  of 
spacecraft  mass.  Propulsion  requirements  for  these  small  satellites  can  vary  significantly  given  the 
diversity  of  small  satellite  missions,  but  are  generally  characterized  by  their  very  low  thrust  values  (- 
mN)  and  low  impulse  bits  for  attitude  control  (~  10"s  p,N-s).  Scaling  traditional  propulsion  systems  down 
in  power  and  size  to  suit  the  needs  of  small  satellites  is  a  major  engineering  challenge. 

Cold  gas  thrusters  are  low-performance  devices  limited  by  low  values  of  specific  impulse  (-50  s).3 
Advantages  of  these  devices  compared  to  other  propulsion  systems  include  their  low  system  complexity, 
low-cost,  and  absence  of  satellite  contamination  problems  when  using  inert  propellants  (e.g.  N2). 
Micronozzle  cold  gas  thrusters  are  suitable  for  attitude  control  applications  of  small  satellites.  They  can 
also  be  used  as  primary  propulsion  system  if  the  Av  requirements  are  less  than  -100  m/s.  Since  nozzle 
throat  diameters  in  the  micrometer  range  are  required  for  micronozzle  propulsion,  the  production  of  these 
microthrusters  rely  on  MEMS  fabrication  technologies.4'6 

The  extremely  small  dimensions  of  microdischarges  combined  with  intense  and  controllable  gas 
heating  can  be  exploited  in  microthrusters  technologies.  Like  conventional  cold  gas  thrusters, 
electrothermal  microdischarge  propulsion  systems  can  have  extremely  low  mass  and  volume  footprints. 
They  are  electrically  simple,  and,  unlike  other  electric  propulsion  technologies,  they  do  not  require 
auxiliary  systems  such  as  neutralizers,  heaters,  or  magnets.  The  thrust  from  these  devices  is  widely 
tunable  by  varying  power  levels.  A  significant  improvement  in  the  specific  impulse  compared  to 
conventional  micronozzle  cold  gas  thrusters  can  potentially  be  achieved  with  these  systems. 

The  Micro  Plasma  Thruster  (MPT)  concept  consists  of  a  direct-current  microdischarge  in  a  geometry 
comprising  a  constant  area  flow  section  followed  by  a  diverging  exit  nozzle.  Geometric  dimensions  are 
on  the  order  of  hundreds  of  microns.  In  the  MPT  concept,  a  stable  microdischarge  is  generated  at  low 
propellant  flow  rates  (few  seem).  The  propellant  gas  in  the  microdischarge  is  eventually  expelled  from 
the  diverging  nozzle.  Thrust  is  produced  by  the  expanding  gas  flow,  which  is  heated  by  the 
microdischarge.  The  MPT  has  potentially  important  advantages  over  competing  concepts  such  as  the 
micro  resistojet  thruster8,  which  also  relies  on  preheating  an  expanding  gas  stream.  The  important 
difference  between  a  resistojet  and  the  MPT  is  the  nature  of  power  deposition  into  the  gas.  In  a  resistojet, 
the  heating  element  is  essentially  the  hottest  part  of  the  system  and  heat  addition  to  the  gas  occurs  via 
conduction/convection  to  the  gas  stream.  There  is  clearly  an  upper  limit  on  the  heating  element 
temperature  which  in  turn  limits  the  amount  of  heat  delivered  to  the  gas.  In  the  MPT,  the  gas  is  heated 
directly  by  the  plasma.  Importantly,  if  the  plasma  thermal  power  deposition  occurs  far  away  from  the 
surfaces,  much  higher  thermal  energies  can  be  delivered  to  the  gas. 

Successful  design  of  the  microthruster  requires  a  detailed  understanding  of  microdischarge  phenomena 
as  well  as  plasma-flow  interactions.  Experimental  investigations  of  the  microdischarge  are  hindered  by 
the  extremely  small  geometries,  where  well-established  plasma  diagnostic  techniques  are  not  applicable. 
High-fidelity  computational  models  can  therefore  have  a  big  impact  in  developing  a  fundamental 
understanding  of  the  MPT. 

In  previous  studies,  Kothnur  and  Raja9  studied  direct-current  microdischarges  in  a  flowing  gas  stream 
for  applications  in  electro-thermal  microthrusters.  For  currents  around  1  mA,  the  microdischarge  is  found 
to  operate  in  an  abnormal  glow  mode  with  positive  differential  resistivity.  An  increase  in  input  electrical 
power  results  in  almost  linear  increase  in  the  gas  temperatures;  this  property  of  microdischarges  is  a  key 
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feature  that  can  be  exploited  in  the  MPT  concept.  Arakoni  et  al 10  studied  a  nominal  microdischarge 
configuration  with  a  flowing  gas  that  resembles  a  thruster,  although  the  flow  modeling  approach  provides 
an  inadequate  description  of  the  gas  dynamics  expansion  through  the  nozzle  and  subsequent  expansion 
into  vacuum.  Gas  temperatures  exceeding  1000  K  are  reported  for  power  densities  of  10’s  kW/cm3  at 
upstream  pressure  of  tens  of  Torr.  They  found  that  the  nozzle  length  and  the  location  of  the  discharge  in 
the  nozzle  have  an  important  influence  on  the  incremental  thrust  (above  that  of  the  cold  flow). 

Section  II  describes  the  mathematical  model  and  numerical  approach  used  in  this  study.  The  model 
comprises  a  plasma  module  coupled  to  a  flow  module.  The  plasma  module  provides  a  self-consistent, 
multi-species,  multi-temperature  description  of  the  microdischarge  phenomena  while  the  flow  module 
provides  a  description  of  the  low  Reynolds  number  compressible  flow  through  the  MPT.  Section  III 
describes  simulation  results  and  a  discussion  of  the  MPT  phenomena.  Finally,  we  conclude  in  Section  IV. 

II.  Methodology 

Coupled  plasma  and  gas  flow  phenomena  are  encountered  in  a  number  of  applications  requiring  the 
development  of  an  integrated  plasma  model  coupled  to  a  gas  dynamic  flow  model.11  In  the  MPT,  the 
microdischarge  plasma  is  in  a  highly  non-equilibrium  state  with  disparate  electron  and  gas  temperatures 
and  non-equilibrium  finite-rate  chemical  effects.  The  gas  flow  is  characterized  by  low-Reynolds  number 
viscous  dominated  effects  combined  with  high  Mach  number  compressibility  effects  owing  to  the 
expansion  of  the  relatively  high-pressure  (-100  Torr)  micronozzle  gas  stream  into  vacuum.  Both  the 
plasma  and  flow  phenomena  are  strongly  coupled:  the  plasma  causes  gas  heating,  which  modifies  the  gas 
density,  and  hence  the  flow  field:  the  gas  flow  velocities,  in  turn,  affects  the  distributions  species  and 
temperature  of  plasma  in  the  discharge.  The  integrated  plasma-flow  model  used  in  this  study  is  described 
below. 


III.  Plasma  Module 


i.  Governing  Equations 

The  densities  of  individual  species  in  the  discharge  are  determined  by  the  species  continuity  equation, 

^-+V.(nkV)+V.fk=Gk  k  =  l,...,Kg(k*kb),  (1) 

where  is  the  number  density  of  species  E  Pis  the  mass-averaged  gas  flow  velocity,  is  the  drift- 
diffusion  flux  of  species  k,  Qk  is  the  gas-phase  species  generation  rate  through  plasma  chemical  reactions, 

and  Kg  is  the  total  number  of  gas  species.  A  dominant  neutral  species  is  identified  as  the  “background” 
species,  and  is  designated  as  fo;  its  density  being  determined  to  satisfy  the  ideal  gas  law, 

p  =  nk„kBTg  +  nekBTe  +  nkkBTg  ,  (2) 

k±ke  ,k±kb 


where  P  is  the  local  pressure,  kB  is  Boltzmann’s  constant,  subscripts  e  and  g  denote  electrons  and  the 
background  gas,  respectively. 

The  electron  temperature,  Te ,  is  explicitly  computed  in  order  to  determine  the  reaction  rates  due  to 
electron  impact  and  the  electron  transport  properties.  For  this  purpose,  the  following  electron  energy 
equation  is  solved. 
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(3) 


f(f  We)  +  V  •  fyBneTeV)  +  V  •  fyBTete)- V  •  (7eVre) 

=  e(fe  +  «/)  •  ^<t>-\kBne  ^-(Te  -  Tg)  Vew  -  e  £  A Ejrj 

1  mh  M 

where  rje  is  the  thermal  conductivity  of  electrons,  me  and  mkh  are  the  molecular  mass  of  electron  and 
dominant  background  gas  respectively,  yek  is  the  electron  momentum  transfer  collision  frequency  with 

the  background  gas,  A Ej  is  the  energy7  lost  per  electron  (in  eV)  in  an  inelastic  collision  event  represented 

by  a  gas-phase  reaction  j,  ry  is  the  rate  of  progress  of  a  reaction  /,  and  Ig  is  the  total  number  of  gas-phase 
reactions.  In  our  model,  we  assume  a  common  temperature  for  all  heavy  species  (ions  and  neutrals),  this 
temperature  being  solved  for  in  the  flow  module. 

The  self-consistent  electrostatic  potential  is  determined  using  the  Poisson’s  equation, 

VV=--ZZ»',»-  (4) 

k 

where  (j)  is  the  potential,  e  is  the  unit  electric  charge,  £q  is  the  permittivity  of  free  space,  and  ZK  is  the 
charge  number  of  species  k  (e.g.,  -1  for  electrons). 

The  ion  flux  is  evaluated  as  rijV  +  f)  =  njUj,  where  is  the  ion  fluid  velocity,  which  is 

determined  using  the  first  moment  of  the  ion  species  Boltzmann  equation,  i.e.,  the  ion  species  momentum 
equation, 


drijUj 

dt 


+  V  •  ( njUjUj )  = 


v  (/,  -  -Mh  vM/  -  m  (uj  -v)Vj. 

m i  mi 


(5) 


The  first  term  on  the  left-hand  side  of  Eq.  (5)  is  the  unsteady  term  and  the  second  term  is  the  ion  inertia 
term  which  is  non-negligible  as  the  pressure  approaches  vacuum  conditions  near  the  MPT  exit.  The  first 
term  on  the  right-hand  side  of  Eq.  (5)  is  the  electrostatic  force  on  the  ion  fluid,  the  second  is  the  gradient 
diffusion  of  ion  fluid  momentum,  and  the  last  term  is  the  loss  of  ion  fluid  momentum  owing  to  collisions 
(friction)  with  other  species.  The  last  term  assumes  that  the  ion  density  is  small  compared  to  the 
dominant  background  species  density,  whose  fluid  velocity  corresponds  to  the  mean-mass  velocity  of  the 
gas.  In  Eq.  (5),  7}  is  the  ion  temperature,  which  is  assumed  equal  the  local  heavy  species  (gas) 

temperature  Tg  and  V\  is  the  ion  momentum-transfer  collision  frequency.  The  mass-averaged  flow 
velocity  and  the  gas  temperature  are  solved  for  using  the  flow  module  discussed  in  the  next  section. 

ii.  Transport  Formulation 

For  species  other  than  the  ions  (neutrals  and  electron),  species  number  flux  is  evaluated  using  the  drift- 
diffusion  approximation,  given  as. 


f  k  =  ~Mknk^ </>  ~  Dk^nk  >  (6) 

where,  -  Z^elmj^v^  is  the  charged  species  mobility  and  D £  ^kgT^/m^v^  is  the  species  diffusion 
coefficient.  The  electron  thermal  conductivity  appearing  the  electron  energy  equation  is  given  as 
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Tje  =  (5/2)  ne  De .  The  operating  (propellant)  gas  considered  for  this  study  is  argon.  The  electron  transport 
properties  are  determined  apriori  by  solving  a  zero-dimensional  electron  Boltzmann  equation 
(“BOLSIG+”12)  with  appropriate  electron  energy-dependent  cross  sections  for  a  range  of  reduced  electric 
fields  (E/N).  The  electron  properties  and  the  electron  mean  energies  are  determined  as  a  function  of  EJN 
by  the  Boltzmann  solver.  These  results  are  then  used  to  tabulate  the  properties  as  a  function  of  the 
electron  mean  energy  (temperature)  for  use  in  the  plasma  discharge  simulations.  The  ion  transport 
properties  are  derived  from  experimental  mobility  data13  based  on  which  the  collision  frequency  vz*  can 

be  determined.  Hard-sphere  cross  sections  phased  on  the  Lennard-Jones  interaction  potentials14  are 
used  to  determine  the  neutral  transport  properties  using  vk  =  ,  where  g ^  is  the  relative  thermal 

speed  of  species  k. 


iii.  Boundary  Conditions 

For  the  electron  continuity  equation,  the  electron  flux  at  solid  walls  is  specified  using  a  kinetic 
Maxwellian  flux  conditions  for  consumption  of  electrons  at  the  walls  with  a  return  flux  of  secondary 
electrons  from  the  walls,  i.e. 


F,  -h  =  —n 


W1/2 

j 


■rrAr+-t’ 


(7) 


where,  h  is  the  unit  normal  vector  pointed  toward  the  wall.  The  solid  walls  include  the  electrode  and 
dielectric  surfaces.  The  first  term  in  Eq.  (7)  is  the  Maxwellian  flux  of  electrons  to  the  surface  and  the 
second  term  is  the  secondary  electron  emission  flux  from  the  surface.  In  this  study,  we  assume  that  the 
secondary  emission  coefficient,  y ,  includes  contributions  from  ion  impact,  fast  atoms,  metastable  atoms, 
dimmers,  and  ultraviolet  photons. 

For  the  ions  continuity  and  ions  momentum  equations,  the  flux  of  positive  ions  is  set  to  0  for  all  solid 
boundaries  for  which  the  electric  field  points  away  from  the  wall  (E-n<  0).  If  the  electric  field  points 
towards  the  wall  (£■«  >0),  the  ion  flux  is  extrapolated  from  the  interior  by  imposing  a  zero  gradient  of 
the  ion  flux  at  the  boundary. 

The  Maxwellian  flux  condition  is  imposed  for  the  neutral  species,  given  by 


kjTg 

nmn  j 


(8) 


For  the  electron  energy  equation,  the  following  energy  flux  is  imposed  at  solid  walls. 


5_ 

2' 


q:  =-kBTere 


w 


(9) 


where,  r Y  ls  the  electron  wall  number  flux. 

The  potential  on  dielectric  surfaces  is  determined  using  Gauss’  law  for  the  dielectric-gas  interface 
boundary  condition  which  depends  on  the  total  accumulated  surface  charge  density  at  the  surface.  9  The 
equation  for  evolution  of  the  net  surface  charge  density,  ps  is  given  by, 
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(10) 


fy>s 

8t 


i  eZkYk-h. 
k= 1 


iv.  Plasma  Chemistry  and  Surface  Processes 

A  pure  argon  plasma  gas  chemistry  is  used  which  comprises  six  species:  electrons  (e),  atomic  argon 
ions  ( Ar+ ),  molecular  argon  ions  ( Ar2+ ),  electronically  excited  atoms  (Ar*),  electronically  excited 
molecules  (Ar2  ),  and  the  background  argon  atoms  (Ar).  Dimer  species  are  included  because  of  the 
relatively  high  pressures  in  the  upstream  section  of  the  micronozzle  (—100  Torr).  The  list  of  reactions 
considered  in  the  study  is  given  in  Table  I  and  comprises  electron  impact  ionization  and  excitation 
reactions.  Penning  ionization  reactions,  three-body  reactions  for  dimer  excited  species  and  ion  formation, 
quenching  and  de-excitation  reactions.  As  in  the  case  of  electron  transport  properties,  the  electron  impact 
reaction  rates  are  determined  apriori  by  solving  the  zero-dimensional  electron  Boltzmann  equation 
(“BOLSterE’12)  with  appropriate  electron  energy-dependent  reaction  cross  sections.  These  rates  are 
tabulated  as  a  function  of  the  electron  temperature  using  the  same  procedure  as  for  the  electron  transport 
properties. 

At  solid  surfaces  all  excited  species  and  charged  species  are  assumed  to  get  quenched  with  unity 
sticking  coefficient.  Upon  quenching  at  surfaces,  each  dimer  ion  and  excited  species  is  assumed  to  return 
to  plasma  as  a  pair  of  ground  state  neutral  Ar  atoms,  while  the  monomer  species  return  as  single  Ar  atoms. 

IV.  Flow  Module 


Conservation  equations  for  the  gas  mass  density,  mass-averaged  gas  velocity,  and  the  gas  energy  are 
solved  using  compressible  Navier-Stokes  equations  in  axisymmetric  form.  These  equations  can  be 
written  as 


+  V  •  Finviscid  —  V  •  Fvjsc  4-  S  ?  (11) 

at 

where,  U  represents  the  conservative  variables  (p,  pVx ,  pVr ,  pet)  (p  is  the  background  density,  Vx  is  the 

axial  velocity,  Vr  is  the  radial  velocity,  and  et  is  the  total  internal  energy  of  the  gas),  Finvi$cid  represents 

the  inviscid  flux  terms,  Fvjsc  represents  the  viscous  flux  terms,  and  S  represents  the  source  terms. 

The  flow  solution  influences  the  plasma  discharge  through  the  pressure,  temperature  and  velocity 
fields  appearing  in  the  plasma  governing  equations  (Eqs.  (l)-(5)).  On  the  other  hand,  electrostatic  forces 
and  electrothermal  heating  act  as  external  body  forces  and  external  heat  source,  respectively,  in  the 
Navier-Stokes  equations.  These  source  terms  are  defined  as 

7  =  21  -eneE  .  (12) 


for  the  electrostatic  force,  and 


Sheat  •  Vfl  +  4 kBn,  -e- 
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(13) 


for  the  electrothermal  heating,  such  that  s  =  (o,  fx,  fr,  Sheatf-  The  first  term  on  the  RHS  of  Eq.  (13) 
represents  the  effective  ion  Joule  heating,  taken  here  to  be  a  fraction  aj  (0  <  aj  <  1 )  of  the  local  ion 
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Joule  heating  (  Jion  •  E ,  where  Jion  is  the  total  ion  current  density  and  E  is  the  local  electric  field).  This 

fraction  is  <  1  for  cases  where  the  ion  mean  free  path  is  comparable  or  larger  than  the  characteristic 
length  scale  for  the  plasma,  which  results  in  incomplete  conversion  of  the  local  kinetic  energy  gained  by 
the  ions  in  the  electric  field  to  ion/heavy  species  thermal  energy.  For  microdischarge  conditions,  we  use  a 
fixed  value  of  aj  =0.25  as  suggested  in  Ref.  15.  The  Joule  heating  terms  appearing  in  the  electron 
energy  equation  (Eq.  (3))  and  in  Eq.  (13)  are  evaluated  using  a  species  flux  reconstruction  approach.16 

V.  Numerical  Approach 

Both  the  plasma  governing  equation  and  the  compressible  flow  Navier-Stokes  equations  are  spatially 
discretized  using  a  cell-centered  finite  volume  approach  on  an  unstructured  mesh  with  mixed  mesh  cell 
types.  A  steady  state  solution  is  sought  in  all  cases.  Both  the  plasma  and  flow  governing  equations  are 
solved  as  transient  problems  with  time-stepping  of  the  solution  to  a  steady  state.  The  plasma  governing 
equations  are  solved  using  an  implicit  time-discretization  approach  with  local  linearization  of  the 
governing  equations  at  each  time  step.  For  the  flow  governing  equations,  the  inviscid  flux  terms  are 
evaluated  with  the  Advection  Upstream  Splitting  Methods  (AUSM)17  and  the  viscous  flux  terms  are 
evaluated  using  the  Haselbacher  approach.18  The  flow  equations  are  also  solved  using  implicit  time- 
discretization,  with  a  dual-time  stepping  approach  to  iteratively  solve  for  the  solution  at  each  time  step. 
The  viscous  terms  in  the  flow  equations  require  computation  of  the  solution  variable  gradients  at  cell 
centers.  We  use  a  Green-Gauss  method  to  reconstruct  these  gradients  based  on  cell-centered  values  of 
these  variables.  Finally,  the  compressible  flow  can  develop  shock  discontinuities.  The 
Venkatakrishnanw  flux  limiter  approach  is  used  to  stabilize  and  produce  monotone  solutions  in  the 
presence  of  such  discontinuities. 


VI.  Physical  Operating  Conditions 

Figure  1  shows  the  geometry  of  the  MPT  used  for  the  simulations.  The  geometry  consists  of  an 
axisymmetric  constant  area  “pipe”  section  of  500  pm  length,  followed  by  a  diverging  section  that  is  200 
pm  in  length,  which  is  terminated  by  a  150  pm  long  constant  area  section.  The  radius  of  the  upstream 
constant  area  pipe  section  is  50  pm,  and  the  exit  section  is  150  pm  in  radius.  The  ring  shaped  electrodes 
have  an  axial  thickness  of  150  pm,  while  the  dielectric  layer  has  an  axial  thickness  of  550  pm.  The  mesh 
consists  of  about  3600  cells,  which  includes  a  combination  of  triangles  and  quadrilaterals. 

The  flow  direction  is  from  left  to  right  in  Fig.  1.  The  flow  enters  the  constant  area  pipe  section  on  the 
left  and  exits  the  domain  along  the  curved  (arc-shaped)  boundary  on  the  right.  Four  boundary  sections  are 
not  modeled  as  solid  walls  for  the  gas  discharge  governing  equations.  Zero-flux  boundary  conditions  are 
imposed  at  the  boundary  section  formed  by  the  symmetry  axis.  For  numerical  stability  reasons,  zero-flux 
boundary  conditions  are  also  used  for  the  gas  discharge  governing  equations  at  the  inlet  and  outer  cathode 
surface  boundaries.  These  boundaries  are  sufficiently  away  from  the  main  discharge  region  and  do  not 
influence  the  results.  The  farfield  boundary  on  the  right  side  of  the  computational  domain  is  modeled  as 
an  “outflow”  boundary,  where  the  plasma  variables  (species  number  densities  and  electron  energy 
density)  are  convected  away  by  the  gas  flow,  while  zero-flux  boundary  conditions  are  imposed  for  the 
Poisson’s  equation.  The  value  of  the  plasma  variable  from  the  adjacent  interior  cell  is  used  to  interpolate 
the  flux  at  this  farfield  boundary.  The  power  input  is  provided  by  applying  a  fixed  positive  DC  voltage  at 
the  anode  (without  ballast  resistance),  while  the  cathode  is  grounded.  For  the  base  case,  the  secondary 
electron  emission  coefficient  is  set  at  a  value  of  0.03,  corresponding  to  a  nickel  surface  interacting  with  an 
argon  plasma22. 

The  inlet  total  (stagnation)  pressure  is  100  Torr  (for  the  base  case),  and  a  small  but  non-zero  outlet 
pressure  (0.05  Torr)  is  required  to  stabilize  the  numerical  scheme  in  the  “vacuum”  part  of  the  domain. 
The  inlet  gas  static  temperature  is  fixed  at  300  K.  The  solid  wall  temperatures  are  fixed  at  300  K.  The 
inlet  flow  velocity  is  computed  self-consistently  using  the  inlet  flow  total  conditions,  the  inlet  static 
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pressure  extrapolated  from  the  interior  of  the  domain  and  isentropic  expansion  relations.  No-slip 
boundary  conditions  are  applied  at  solid  walls  for  the  momentum  equation.  The  Knudsen  numbers  for  the 
flow  are  sufficiently  low  over  much  of  the  MPT  device  length  (about  0.01  at  the  inlet  to  about  0.08  at  the 
exit  plane)  that  slip-flow  and  wall-temperature  jump  boundary  condition  is  not  necessary.  This  has  been 
verified  by  running  simulations  with  the  appropriate  jump  boundary  conditions. 

VII.  Results  and  Discussion 


A.  Base  Case 

Table  II  lists  discharge  conditions  for  which  results  are  presented.  For  the  base  case,  the  discharge 
voltage  is  set  at  a  value  of  750  V,  and  the  computed  discharge  current  is  0.87  mA.  The  total  stagnation 
pressure  (100  Torr)  and  static  temperature  (300  K)  imposed  at  the  inlet  determines  inlet  velocity  which  is 
computed  at  about  100  m/s  corresponding  to  a  flow  rate  of  5.2  seem.  Plasma  properties  at  steady  state  are 
shown  in  Fig.  2  for  the  base  case.  The  space  charge  density  is  high  enough  that  the  bulk  plasma  extends 
well  into  the  diverging  part  of  the  micronozzle  creating  a  hollow-cathode-like  annular  cathode  sheath. 
Under  these  operating  conditions,  the  cathode  sheath  is  about  100-pm  thick  (see  Fig.  2(a)),  occupying  a 
significant  fraction  of  the  total  discharge  volume.  In  the  cathode  fall,  the  potential  drops  by  nearly  750  V 
over  -100  |im  producing  a  characteristic  field  strength  of  -75  kV  cm'1  and  a  reduced  electric  field 
(electric  field  /  gas  number  density)  of  ~105  Td  (Townsend).  Electron  density  contours  are  presented  in 
Fig.  2(b)  and  show  two  local  maxima.  The  first  peak  is  inside  the  constant  area  pipe  section  (-7x1 0W 
#/m3),  and  the  second  peak  is  in  the  diverging  section  of  the  nozzle  (-3x10°°  #/m3).  The  well-defined 
cathode  sheath  structure  observed  in  Fig.  2(a)  is  apparent  in  the  density  contours  with  an  abrupt  drop  in 
the  electron  densities  in  the  diverging  micronozzle  section.  The  MPT  discharge  clearly  operates  in  the 
glow  discharge  mode  rather  than  the  Townsend/predischarge  mode.  The  electron  temperature  contour  is 
shown  in  Fig.  2  (c).  The  electron  temperature  is  only  shown  in  regions  of  the  discharge  where  the 
electron  number  density  is  greater  than  3x1 01  #/nf,  i.e.  10"3  of  the  peak  value  of  the  electron  number 
density.  The  electron  energy  content  is  negligible  in  the  rest  of  the  domain,  although  the  electron 
temperatures  can  increase  exponentially  in  these  regions  of  vanishingly  small  electron  densities  owing  to 
a  numerical  artifact  of  the  fluid  model.  The  electron  temperature  remains  nearly  uniform  («2  eV)  in  the 
bulk  plasma  over  most  of  the  constant  area  pipe  section,  and  gradually  increases  in  the  diverging  section 
of  the  devices  to  temperatures  around  3  eV  at  the  cathode  sheath  edge. 

The  generation  rate  of  electrons  in  the  discharge  through  gas-phase  reactions  is  shown  in  Fig.  2(d). 
Significant  generation  of  electrons  is  observed  over  the  entire  micronozzle  with  maximum  generation 
observed  in  the  diverging  section.  This  location  corresponds  to  a  relatively  higher  electron  temperature 
(^3  eV)  and  reasonably  high  background  densities  in  the  expanding  gas  flow.  The  argon  monomer  ion 
(Ar+)  and  dimer  ion  (Ar2+)  number  densities  are  shown  in  Figs.  2(e)  and  2(f),  respectively.  Atomic  argon 
ions  constitute  the  dominant  ion  species  in  the  microdischarge.  High  background  densities  favor  the 
three-body  reaction  Gj2  (see  Table  I)  that  forms  dimer  species.  Therefore,  most  dimer  ions  are  located  in 
the  constant  area  pipe  section,  where  the  pressure  is  relatively  high.  The  dielectric  surfaces  support  a  net 
negative  charge  owing  to  electron  trapping  which  in  turn  supports  a  positive  sheath  over  the  entire 
dielectric  length  between  the  anode  and  the  cathode.  This  is  apparent  from  the  excess  of  positive  ions 
compared  to  electrons,  just  above  the  dielectric  surfaces.  The  net  negative  charge  at  the  dielectric 
enforces  the  electron  wall  flux  to  equal  the  total  positive  ion  wall  flux  (at  steady  state),  so  no  net  current  is 
drawn  through  the  dielectric  surfaces. 

The  Ar  ion  number  flux,  computed  using  the  ion  momentum  equation  (Eq.  (5))  is  shown  in  Fig.  3. 
Most  of  the  ions  are  generated  in  a  localized  region  close  to  the  axis  in  the  diverging  portion  of  the  device 
and  subsequently  drift  and  diffuse  away  from  this  region.  Most  of  the  ions  eventually  drift  towards  the 
cathode  and  the  dielectric  surface  in  the  vicinity  of  the  cathode,  where  they  are  quenched.  A  negligibly 
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small  fraction  of  the  ions  (<  1%)  leave  the  nozzle.  Consequently,  the  electrostatic  ion  thrust  contribution 
to  the  overall  thrust  produced  by  the  MPT  device  is  negligible. 

The  base  case  simulation  was  also  run  with  the  drift-diffusion  approximation  (Eq.  (6))  used  for  ion 
transport,  instead  of  the  ion  momentum  equation  (Eq.  (5)).  Comparison  of  the  discharge  results  for  the 
drift-diffusion  ion  simulation  case  (not  shown)  with  the  ion  momentum  case  indicate  that  the  results  can 
be  significantly  different  (e.g.  the  peak  value  of  the  monomer  ion  number  densities  differs  by  about  20% 
between  the  two  cases),  confirming  the  non-negligible  role  of  the  ion  inertia  in  the  MPT  configuration. 

The  flow  field  solution  for  the  base  case  is  shown  in  Fig.  4.  The  major  influence  of  the  plasma  on  the 
flow  field  is  through  the  heat  addition  that  is  shown  in  Fig.  4(a).  We  have  confirmed  through  additional 
simulations  that  the  electrostatic  forcing  source  term  in  the  flow  momentum  equation  (Eq.  (12))  has  a 
negligible  effect  on  the  overall  flow  and  discharge  solutions.  The  thermal  power  deposition  to  gas  (Eq. 
(13))  originates  in  large  part  from  ion  Joule  heating  in  the  cathode  sheath.  A  peak  in  the  gas  heating  is 
observed  close  to  the  MPT  surface  at  the  comer  between  the  dielectric  and  the  cathode  where  the  electric 
fields  are  the  highest.  The  gas  temperature  contours  are  shown  in  Fig.  4(b).  The  gas  temperature 
contours  show  a  peak  value  of  about  950  K.  The  location  of  the  peak  temperature  is  slightly  downstream 
and  closer  to  the  axis  compared  to  the  location  of  peak  thermal  power  deposition.  The  localized  power 
deposition  in  the  immediate  vicinity  of  the  cathode  surface  leads  to  significant  heat  loss  to  the  solid  wall, 
indicating  that  thruster  performance  could  be  improved  in  the  future  by  optimizing  the  design  to  minimize 
this  loss. 

The  axial  and  radial  velocity  components,  Mach  number,  and  pressure  fields  are  shown  in  Figs.  4(c)  to 
4(f).  The  axial  velocity  increases  rapidly  at  the  start  of  the  diverging  section  owing  to  gas  dynamic 
expansion,  followed  by  a  small  decrease  owing  to  compression  waves  launched  from  the  straight  cathode 
section.  Subsequently,  the  axial  velocity  increases  once  it  expands  out  the  exit  plane.  Peak  axial 
velocities  greater  than  550  m/s  are  seen  outside  the  exit  plane.  The  Mach  number  contours  follow  the 
trends  in  the  axial  velocity,  with  an  increase  at  the  start  of  the  diverging  section,  followed  by  a  decrease 
near  the  cathode  section  (owing  to  axial  velocity  decrease  as  well  as  an  increase  in  the  sound  speed  due  to 
high  temperatures),  and  a  subsequent  increase  beyond  the  exit  plane.  The  radial  velocity  profiles  reflect 
the  mass  continuity  requirements  for  the  gas  expansion  process.  The  pressure  profile  shows  a  gradual, 
almost  linear  Poiseuille,  pressure  drop  in  the  straight  pipe  section  of  the  micronozzle,  followed  by  a  rapid 
decrease  in  the  diverging  section  and  out  the  exit  plane.  We  must  note  that  no  sharp  shock  discontinuities 
are  observed  in  the  gas  flow  solution,  owing  mainly  to  the  very  low  Reynolds  number  viscous  dominated 
nature  of  the  flow,  despite  the  high  Mach  numbers.  For  these  base  case  operating  conditions,  the 
computed  thrust  is  100  pN  (compared  to  67  pN  for  the  cold  gas  micro  nozzle),  which  corresponds  to  a 
specific  impulse  of  74  s.  The  electrostatic  component  of  thrust  (owing  to  ion  momentum  leaving  the  exit 
plane)  is  found  to  be  negligible  and  the  thrust  is  almost  entirely  due  to  neutral  momentum  transport  out 
the  exit  plane. 

Figure  5  shows  a  schematic  of  the  energy  flow  pathways  in  the  MPT.  The  electrical  power  input  (the 
discharge  current  times  the  discharge  voltage)  for  the  base  case  is  650  mW.  Of  this,  most  of  the  power 

(about  550  mW)  is  attributed  to  ion  Joule  heating  integrated  over  the  entire  domain  (  J  Jlon  • EdV ),  while 

the  electron  Joule  heating  ( J  • EdV )  is  about  45  mW.  The  mismatch  between  the  input  electrical 

power  and  the  integrated  Joule  heating  in  the  domain  (ion  and  electron)  is  attributed  to  numerical 
discretization  error.  According  to  our  model  a  large  fraction  of  the  ion  Joule  heating  is  lost  to  the  wall 
while  the  remaining  is  thermalized  to  provide  a  source  term  to  the  gas  energy  equation.  This  fraction  is 
set  by  the  parameter  aj  which  equals  0.25  in  these  simulations.  Consequently,  75%  of  the  intrinsic  ion 

Joule  heating  is  lost  directly  to  the  wall  surfaces.  A  small  fraction  (less  than  1%)  of  the  power  deposited 
in  the  heavy  species  thermal  pool  (gas  energy)  is  lost  through  inelastic  collisions,  while  the  rest  is  lost  at 
solid  walls  and  at  the  outflow  section.  The  rate  of  energy  transfer  from  the  electron  thermal  pool  to  the 
heavy  species  (gas)  thermal  pool  through  inelastic  collisions  amounts  to  a  small  fraction  of  the  electron 
Joule  heating  (less  than  1%).  A  significant  amount  of  power  (35  mW  or  ~80%  of  the  electron  Joule 
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heating)  is  lost  through  inelastic  collisions,  while  the  rest  of  the  power  deposited  in  the  electron  thermal 
pool  is  lost  through  transport  at  solid  walls  and  at  the  outflow  section 

VIII.  Effect  of  Power  Input 

Figure  6  shows  the  effects  of  power  input  on  discharge  characteristics.  Results  are  presented  for  the 
base  case  (where  the  discharge  voltage  is  fixed  at  750  V)  and  for  a  case  where  the  discharge  voltage  is 
fixed  at  1000  V.  For  both  cases,  the  inlet  total  pressure  is  100  Torr  and  the  flow  rate  is  5.2  seem.  The 
plasma  is  more  intense  (i.e.  the  ionized  species  number  densities  are  higher)  for  higher  power  inputs. 
Electron  and  dominant  ion  Ar  number  densities  increase  by  a  factor  of  ~2  as  the  power  input  is  increased 
from  650  mW  (base  case)  to  1800  mW. 

Figure  7  shows  the  effect  of  power  input  on  flow  properties.  Contours  of  the  power  deposition  in  the 
MPT  are  shown  in  Fig.  7(i).  The  net  power  deposition  into  neutral  gas  scales  almost  linearly  with  the 
power  input  as  the  ratio  between  these  two  quantities  remains  equal  to  about  20%  for  both  cases.  Gas 
temperatures  contours  are  shown  in  Fig.  7(ii).  A  peak  temperature  of  1850  K  is  reached  for  the  high 
power  case  (compared  to  950  K  for  the  base  case).  Importantly,  these  results  show  that  changing  the 
external  power  input  is  an  effective  method  to  control  the  level  of  gas  heating  in  the  discharge,  which,  in 
turn,  affects  the  gas  temperature  and  the  thrust  produced  by  the  device.  This  constitutes  a  key  feature  of 
our  proposed  MPT  concept.  The  computed  thrust  increases  to  a  value  of  128  pN  for  the  high  power  case, 
which  corresponds  to  a  specific  impulse  of  95  s. 

IX.  Effect  of  Pressure  /  Flow  Rate 

We  now  vary  the  inlet  total  (stagnation)  pressure  to  200  Torr  (the  base  case  inlet  total  pressure  is  100 
Torr).  Other  conditions  such  as  the  inlet  static  temperature  (300  K)  and  the  discharge  voltage  drop  (750 
V)  are  kept  the  same  as  the  base  case.  The  increase  of  the  inlet  total  pressure  causes  an  increase  of  flow 
rate  through  the  micronozzle  to  a  value  of  1 3  seem  (compared  to  5.2  seem  for  the  base  case).  Contours  of 
the  axial  velocity  are  shown  in  Fig.  8(i).  Higher  flow  velocities  and  gas  mass-densities  (not  shown)  are 
observed  throughout  the  nozzle  for  the  high  pressure  case.  The  change  in  flow  field  has  a  profound  effect 
on  the  microdischarge  characteristics.  The  high  velocities  combined  with  decreased  diffusive  transport 
properties  lead  to  a  further  downstream  convection  of  the  ionized  species.  The  discharge  also  appears 
more  constricted  for  the  high  pressure  case  with  a  slightly  lower  peak  charge  density  for  the  high  pressure 
case.  The  higher  pressures  decrease  the  plasma  conductivity  resulting  in  a  decrease  in  the  discharge 
current  to  0.75  mA  (compared  to  0.87  mA  for  the  base  case).  The  resulting  input  power  decreases  for  the 
high  pressure  case  to  about  560  mW  (compared  to  650  mW  for  the  base  case).  The  peak  gas  temperature 
also  decreases  compared  to  the  base  case  (600  K  compared  to  950  K,  not  shown)  owing  to  lower  power 
deposition  and  also  because  of  the  increased  heat  capacity  of  the  higher  pressure  gas.  Case  studies  were 
also  performed  with  a  lower  total  inlet  pressure  of  60  Torr.  For  these  conditions,  the  microdischarge 
could  not  be  sustained  for  discharge  voltages  ranging  from  750  V  to  1000  V. 

X.  Sensitivity  to  Model  Parameters 


1.  Cathode  Temperature 

Since  a  large  fraction  of  the  input  power  is  lost  to  the  walls,  the  cathode  temperature  can  be  expected 
to  be  significantly  higher  than  the  constant  300  K  assumed  in  the  above  results.  An  accurate  estimate  of 
the  cathode  temperature  would  require  solving  an  energy  balance  equation  for  the  entire  MPT  device 
(including  solid  materials). 

The  effect  of  the  cathode  temperature  on  the  MPT  discharge  parameters  and  performance  is  studied  by 
comparing  the  base  case  solution  to  a  case  where  the  cathode  temperature  is  fixed  at  1000  K.  All  other 
conditions  are  kept  the  same  as  the  base  case.  Figure  9(i)  plots  the  gas  temperature  contours  for  the  two 
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cases.  The  spatial  distribution  of  gas  temperature  appears  strongly  dependent  on  the  cathode  temperature. 
For  the  case  where  the  cathode  temperature  is  fixed  at  1000  K,  the  peak  gas  temperature  reaches  1 150  K 
(compared  to  950  K  for  the  base  case)  and  the  region  of  high  gas  temperatures  occupies  a  larger  fraction 
of  the  MPT  volume  than  for  the  base  case.  Increasing  the  cathode  temperature  reduces  the  temperature 
gradients  at  the  cathode  surface,  thereby  decreasing  the  net  thermal  loss. 

The  axial  velocity  contours  are  shown  for  the  two  cases  in  Fig.  9(ii).  As  the  cathode  temperature  is 
increased,  the  reduced  heat  losses  allow  for  an  increased  expansion  of  the  gas  in  the  diverging  section  of 
the  micronozzle.  The  computed  thrust  therefore  increases  from  a  value  of  100  pN  for  the  base  case  to  a 
value  of  1 12  pN  for  the  1000  K  cathode  temperature  case.  It  should  be  noted  that,  for  this  higher  cathode 
temperature  case,  the  power  input  has  decreased  to  a  value  of  540  mW,  compared  to  650  mW  for  the  base 
case.  The  higher  gas  temperatures  (and  therefore  lower  background  Ar  densities)  result  in  a  decrease  in 
ionization.  The  lower  power/current  can  be  attributed  to  the  lower  charge  densities  and  hence  to  the 
lower  plasma  conductivity. 

2.  Fraction  of  Ion  Joule  Energy  Thermalized  with  the  Gas 

For  high-pressure  microdischarges  and  conventional  discharges  at  very  low  pressures,  a  large  fraction 
of  the  kinetic  energy  of  the  ions  is  deposited  directly  onto  discharge  bounding  surfaces,  particularly  the 
cathode.  This  phenomenon  has  been  simulated  by  Revel  et  al 21  using  a  particle-in-cell  Monte-Carlo 
approach  for  a  one-dimensional  argon  discharge.  For  current  densities  of  ~20  A/m2,  an  interelectrode 
distance  of  1.5  cm,  and  a  background  pressure  of  1  Torr,  they  observe  that  -25%  of  the  ion  Joule  energy 
is  converted  to  gas  thermal  energy,  with  the  remainder  being  transferred  directly  to  the  cathode.  Boeuf  et 
a/15  used  a  nominal  fixed  value  of  aj~  0.25  for  simulation  studies  of  a  100  Torr  xenon  micro-hollow 

cathode  discharge,  which  is  the  same  value  used  in  our  study. 

The  sensitivity  of  our  model  results  to  this  parameter  is  determined  by  varying  aj  to  0.5.  The  power 

deposition  and  gas  temperatures  contours  are  shown  in  Fig.  10  for  both  cases.  Since  most  of  the  thermal 
(gas)  power  depositon  in  the  MPT  is  attributed  to  ion  Joule  heating,  the  thermal  power  deposition 
increases  significantly  from  a  value  of  140  mW  for  the  base  case  to  260  mW  when  aj  =  0.5.  At  the 

same  time,  the  total  input  power  to  the  MPT  decreased  to  540  mW  (compared  to  650  mW  for  the  base 
case).  As  described  in  the  previous  section,  the  increase  of  the  gas  temperature  (to  a  peak  value  of  1350 
K)  results  in  decreased  neutral  density,  with  a  resulting  decrease  in  ionization  rate  and  hence  the  charge 
density.  Consequently,  the  plasma  conductivity  decreases  resulting  in  decreased  current  and  total  power 
deposition  into  the  discharge.  The  increased  fraction  of  power  deposited  into  gas  improves  the  efficiency 
of  the  device.  The  computed  thrust  increases  to  1 16  pN  (compared  to  100  pN  for  the  base  case). 

Our  model  assumes  a  constant  value  of  ay.  In  reality,  this  parameter  will  depend  on  the  local 

background  pressure.  Indeed,  as  the  mean  free  path  of  ions  decreases  with  increasing  pressure,  ions  will 
deposit  a  larger  fraction  of  their  kinetic  energy  to  the  background  gas.  In  future  studies  we  will  explore 
more  accurate  strategies  (e.g.  hybrid  models)  to  model  this  effect  consistently. 

3.  Secondary  Electron  Emission  Coefficient 

The  value  of  the  secondary  electron  emission  coefficient,  y ,  constitutes  an  important  uncertain 
parameter  of  the  plasma  model.  This  coefficient  is  sensitive  to  surface  conditions  and  its  value  can  vary 
by  several  orders  of  magnitude  depending  on  the  surface  material  and  the  cleanliness  of  the  surface.  In 
the  above  simulation  results,  we  chose  a  constant  value  of  y=  0.03  (corresponding  to  a  nickel  surface 
interacting  with  an  argon  plasma).22  The  sensitivity  of  our  results  to  this  parameter  is  analyzed  by 
comparing  our  base  case  to  a  case  with  y  =0.01.  The  discharge  current  is  significantly  reduced  for  the 
model  with  the  lower  secondary  electron  coefficient  (to  a  value  of  0.33  mA).  Lower  y  results  in  fewer 
electrons  released  from  the  surface  through  Auger  processes  per  ion  impact,  which  in  turn  results  in  a 
much  weaker  discharge.  This  behavior  is  characteristic  of  all  gamma-mode  discharges  such  as  the  direct- 
current  microdischarge  presented  in  this  study.2"*'24 
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Figure  11  shows  the  ion  number  density  and  gas  temperature  contours  obtained  with  the  two  models 
for  a  discharge  voltage  of  750  V.  While  the  peak  ion  number  density  has  about  the  same  value  (~3><1020 
#/m3)  for  both  cases,  the  discharge  activity  appears  more  constricted  for  the  lower  y  case.  Since  the  ion 
Joule  heating  scales  with  the  current,  a  significant  decrease  in  power  deposition  into  the  neutral  gas  is 
observed  (not  shown).  This  results  in  a  decrease  of  the  gas  temperature  (shown  in  Fig.  ll(ii))  and  a 
decrease  in  the  computed  thrust  to  a  value  of  87  jiN,  compared  to  100  pN  for  the  base  case. 

XI.  Conclusion 

Microdischarge  and  flow  interaction  phenomena  for  micropropulsion  applications  have  been  studied 
using  a  detailed  self-consistent  computational  model.  The  model  consists  of  a  plasma  module  coupled  to  a 
flow  module  and  is  solved  on  a  hybrid  unstructured  mesh  framework.  The  plasma  module  provides  a  self- 
consistent,  multi-species,  multi-temperature  description  of  the  microdischarge  phenomena  while  the  flow 
module  provides  a  description  of  the  low  Reynolds  number  compressible  flow  through  the 
micropropulsion  system. 

The  Micro  Plasma  Thruster  (MPT)  concept  consists  of  a  direct-current  microdischarge  in  a  geometry 
comprising  a  constant  area  flow  section  followed  by  a  diverging  exit  nozzle.  The  coupled  plasma-flow 
simulations  showed  that  the  electrostatic  component  of  thrust  is  negligible  for  our  current  MPT 
configuration.  Electrothermal  heating  is  due  primarily  to  the  ion  Joule  heating  occurring  near  the  cathode 
surfaces.  A  large  part  of  the  input  power  is  therefore  deposited  into  the  walls  rather  than  going  to  neutral 
gas  heating.  Future  work  must  explore  strategies  that  will  allow  for  power  deposition  away  from  the  wall, 
near  the  centerline.  For  example,  radio-frequency  operation  of  the  microdischarge  may  help  reduce  the 
thermal  load  on  the  microdischarge,  and  increase  the  direct  power  deposition  into  the  neutral  gas. 

For  a  discharge  voltage  of  750  V,  a  power  input  of  650  mW,  and  an  argon  mass  flow  rate  of  5.2  seem, 
the  specific  impulse  of  the  device  is  74  s,  a  factor  of  —1.5  increase  compared  to  the  cold  gas  micronozzle. 
For  these  conditions,  charged  species  densities  on  the  order  of  5xl020  m“3  and  peak  gas  temperatures  of 
-1000  K  are  predicted.  The  microdischarge  remains  mostly  confined  inside  the  micronozzle  and  operates 
in  an  abnormal  regime.  Additional  simulation  studies  on  the  MPT  indicate  that  the  power  input  has  a 
strong  influence  on  overall  discharge  properties.  The  net  power  deposition  into  the  neutral  gas  scales  with 
the  power  input,  providing  a  method  of  controlling  the  gas  temperature  and  the  thrust  level  of  the  MPT. 
This  feature  constitutes  an  important  advantage  over  traditional  cold  gas  thrusters. 

The  sensitivity  of  above  results  to  important  uncertain  parameters  of  the  model  has  been  studied. 
These  parameters  are  the  cathode  temperature,  the  fraction  of  ion  Joule  heating  that  is  locally  converted  to 
gas  thermal  energy,  and  the  secondary  electron  coefficient.  A  higher  cathode  temperature  is  found  to 
reduce  the  heat  losses  at  solid  walls  and  is  beneficial  for  the  thruster  performance,  though  there  is 
obviously  an  upper  limit  to  avoid  damage  to  the  cathode.  The  cathode  temperature  could  be  evaluated  by 
solving  the  solid  thermal  conduction  equation  in  the  dielectric  and  electrode  material  to  quantify  the 
overall  thermal  field  in  the  MPT  device.  The  thermal  loads  from  the  plasma  at  the  wall  surfaces  will  drive 
heat  transfer  into  the  device  structure  potentially  resulting  in  “hot  spots”  that  can  be  investigated  through 
the  simulations.  These  coupled  multi-physics  simulation  results  would  in  turn  provide  an  informed 
approach  to  the  better  designs  of  the  MPT  device  from  a  thermal  viewpoint.  Finally,  the  voltage-current 
characteristics  of  the  device  is  found  to  be  very  sensitive  to  the  value  of  the  secondary  electron  coefficient 
used  in  the  model.  Future  studies  of  the  MPT  should  include  a  calibration  of  these  uncertain  parameters 
and  a  validation  of  the  model  with  experimental  results. 
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# 

Reactions 

Reaction  ratea 

Ref. 

G1 

e  +  Ar  — ►  e  +  Ar 

b 

12 

G2 

e  +  Ar  — ►  e  +  Ar* 

b 

12 

G3 

e  +  Ar  — ►  2e  +  Ar+ 

b 

12 

G4 

e  +  Ar*  — >  2e  +  Ar+ 

b 

12 

G5 

e  +  Ar*  — ►  e  +  Ar 

b 

12 

G6 

e  +  Ar+  — ►  Ar* 

4.0  x  10'13  TV0-5 

20 

G7 

2e  +  Ar+  — >  Ar*  +  e 

5.0  x  10‘27  TV4'7  cm6s"1 

20 

G8 

e  +  Ar2+  — * ►  Ar*  +  Ar 

5.38  xlO'8  T;066 

20 

G9 

2Ar*  — ►  Ar+  +  Ar  +  e 

5x  10'10 

20 

G10 

2Ar2*  — ►  Ar2+  +  2  Ar  +  e 

5x  10'10 

20 

Gil 

Ar*  +  2Ar  — >  Ar2*  +  Ar 

1.14  x  10-32  cm6  s“ 

20 

G12 

Ar+  +  2Ar  -►  Ar2+  +  Ar 

2.5  x  10'31  cm6  s’1 

20 

G13 

Ar2*  -►  2Ar 

6.0  x  107  s'1 

20 

G14 

e  +  Ar2*  — ►  2e  +  Ar2+ 

9  x  10‘8  Tc07  exp(-3.66/Te) 

20 

G15 

e  +  Ar2*  — ►  e  +  2Ar 

10‘7 

20 

a  Rate  coefficients  have  units  of  cm  ’  s'1,  unless  mentioned  otherwise.  Electron  temperatures, 

are  in  eV. 

b  Tabulated  rate  coefficient  as  a  function  of  mean  electron  temperature  was  obtained  by  the 
Boltzmann  equation  solver  “BOLSIG+”.12 

Table  I.  High-pressure  argon  plasma  gas-phase  chemistry  used  in  this  study. 
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Case 

Inlet  total  pressure 
(Ton) 

P  D  (Torr 
cm) 

Flow  rate 
(seem) 

Voltage 

(n 

Current 

(mA) 

Power 

(mil) 

Base  case 

100 

5.5 

5.2 

750 

0.87 

650 

Larger  power  input 

100 

5.5 

5.2 

1000 

1.8 

1800 

Larger  pressure  / 
flow  rate 

200 

11.0 

13 

750 

0.75 

560 

Larger  cathode 
temperature 

100 

5.5 

5.2 

750 

0.72 

540 

Larger  aj 

100 

5.5 

5.2 

750 

0.72 

540 

Lower  yeff 

100 

5.5 

5.2 

750 

0.33 

250 

Table  II.  Operating  conditions  considered  for  this  study.  The  characteristic  discharge  dimension  D  is  taken  to 
be  the  thickness  of  the  dielectric  layer  located  between  the  electrodes  (550  pm)  and  the  characteristic  pressure 
is  taken  to  be  the  inlet  total  pressure. 


150  am  550  150/zm  Symmetry  BC  (axisymmetric) 


Figure  1.  Schematic  of  the  MPT  device  and  computational  mesh.  The  geometry  is  cylindrically  symmetric. 
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(a)  Electrostatic  potential  (V) 


(b)  Electron  number  density  (  #  /  m3  ) 


(d)  Electron  generation  rate  (  #  /  m3-s)  ) 


3.064-16  9.56*16  3.0€*17  9.<?*17  3.C6*1S 


(f)  Ar2+  number  density  (  #  /  m3  ) 


*.#  >.?  7*  3 


(c)  Electron  temperature  (eV) 


2.06*T8  6.36*18  2CC*79  6  36- 


(e)  Ar+  number  density  (  #  /  m3  ) 


Figure  2.  Plasma  properties  in  the  MPT.  The  inlet  total  pressure  is  100  Torr  (13.3  kPa)  and  the  mass  flow 
rate  is  5.2  seem.  The  applied  potential  difference  between  the  electrodes  is  750  V  and  the  power  input  is  650 


mW. 


?t-JJ  ■  it  *21 


Figure  3.  Monomer  ( Ar+ )  ion  number  flux  (  #  /  m2-s)  )  in  the  MPT.  The  contours  indicate  the  magnitude  of 
the  ion  number  flux.  The  vectors  indicate  the  direction  of  ion  flux,  which  the  length  of  arrows  being 
proportional  to  the  flux  magnitude.  The  operating  conditions  correspond  to  the  base  case  indicated  in  Fig.  2 

caption. 
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) 


(b)  Gas  temperature  (  K) 


(c)  Axial  velocity  (  m/s  ) 


(d)  Radial  velocity  (  m/s) 


Figure  4.  Flow  properties  in  the  MPT.  The  operating  conditions  are  indicated  in  Fig.  2  caption. 


Figure  5.  Energy  flow  pathways  in  the  MPT. 


3.5 


(i)  Electron  number  density  (  #  /  m3  ) 


(ii)  Ar+  number  density  (  #  /  ra3  ) 

(a)  Power  input:  650  mW  (b)  Power  input:  1800  mW 

Figure  6.  Electron  and  dominant  ion  ( Ar ,+)  number  density  contours  for  different  values  of  the  power  input. 
The  inlet  total  pressure  is  13.3  kPa  (100  Torr)  and  the  flow  rate  is  5.2  seem  for  both  cases. 


(i)  Power  deposition  (  W  /  m3  ) 


(ii)  Gas  temperature  (K) 


(a)  Power  input:  650  mW  (b)  Power  input:  1800  mW 

Figure  7.  Power  deposition  and  gas  temperature  contours  for  different  values  of  the  power  input.  The  inlet 
total  pressure  is  13.3  kPa  (100  Torr)  and  the  flow  rate  is  5.2  seem  for  both  cases. 
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(i)  Axial  velocity  (  m/s  ) 


(ii)  Electron  number  density  (  #  /  m3  ) 

(a)  Total  inlet  pressure:  100  Torr  (b)  Total  inlet  pressure:  200  Torr 

Figure  8.  Axial  velocity  and  electron  number  density  contours  for  different  values  of  the  inlet  total  pressure  / 
flow  rate.  Other  conditions  are  same  as  the  base  case.  The  discharge  voltage  is  fixed  at  750  V  for  both  cases. 


(i)  Gas  temperature  (K) 


(ii)  Axial  velocity  (  m/s  ) 

(a)  Cathode  temperature:  300  K  (b)  Cathode  temperature:  1000  K 

Figure  9.  Gas  temperature  and  axial  velocity  contours  for  different  values  of  the  imposed  temperature  at  the 
cathode.  The  flow  rate  is  5.2  seem  and  the  applied  potential  difference  between  the  electrodes  is  750  V. 
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(ii)  Gas  temperature  (K) 

(a)  aj  =  0.25  (b)  aj  =  0.5 

Figure  10.  Power  deposition  and  gas  temperature  contours  for  different  fractions  of  ion  Joule  energy 
converted  to  gas  thermal  energy.  The  flow  rate  is  5.2  seem  and  the  applied  potential  difference  between  the 

electrodes  is  750  V. 


(ii)  Gas  temperature  (K) 


(a)  y  =  0.03 


(b)  y  =  0.01 


Figure  11.  Dominant  ion  ( Ar ,4)  number  density  and  gas  temperature  contours  for  different  values  of  the 
secondary  electron  emission  coefficient.  The  flow  rate  is  5.2  seem  and  the  applied  potential  difference 

between  the  electrodes  is  750  V. 
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